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ABSTRACT 

Near-IR imaging of the AGB star IRC+10011 (= CIT3) reveals the presence of a bipolar 
structure within the central ~ 0.1 arcsec of a spherical dusty wind. We show that the image 
asymmetries originate from ~ 10~ 4 M© of swept-up wind material in an elongated cocoon 
whose expansion is driven by bipolar jets. We perform detailed 2D radiative transfer calcu- 
lations with the cocoon modeled as two cones extending to ~ 1,100 AU within an opening 
angle of ~ 30°, imbedded in a wind with the standard r~ 2 density profile. The cocoon ex- 
pansion started < 200 years ago, while the total lifetime of the circumstellar shell is ~ 5,500 
years. Similar bipolar expansion, at various stages of evolution, has been recently observed in 
a number of other AGB stars, culminating in jet breakout from the confining spherical wind. 
The bipolar outflow is triggered at a late stage in the evolution of AGB winds, and IRC+ 1 00 1 1 
provides its earliest example thus far. These new developments enable us to identify the first 
instance of symmetry breaking in the evolution from AGB to planetary nebula. 

Key words: circumstellar matter — dust — infrared: stars — radiative transfer — stars: 
imaging — stars: individual: IRC+1001 1 — stars: AGB and post-AGB 



1 INTRODUCTION 

The transition from spherically symmetric Asymptotic Giant 
Branch (AGB) winds to non-spherical Planetary Nebulae (PNe) 
represents one of the most intriguing problems of stellar astro- 
physics. While most PNe show distinct deviations from spherical 
symmetry, their progenitors, the AGB stars, are conspicuous for 
the sphericity of their winds (see, e.g., review by Balick & Frank 
2002). There have been suggestions, though, that deviations from 
sphericity may exist in some AGB winds, and perhaps could be 
even prevalent (Plez & Lambert 1994, Kahane at al. 1997). Thanks 
to progress in high resolution imaging, evidence of asymmetry has 
become more conclusive for several objects in recent years (V Hya: 
Plez & Lambert 1994, Sahai et al. 2003a; X Her: Kahane & Jura 
1996; IRC+10216: Weigelt et al. 1998 & 2002, Haniff & Buscher 
1998, Skinner et al. 1998, Osterbart et al. 2000; RV Boo: Bergman 
et al. 2000, Biller et al. 2003; CIT6: Schmidt et al. 2002). 

The star IRC+10011 (= IRAS 01037+1219, also known as 
CIT3 and WXPsc), an oxygen-rich long-period variable with a 
mean infrared variability period of 660 days (Le Bertre 1993), is 
one of the most extreme infrared AGB objects. This source served 
as the prototype for the first detailed models of AGB winds by 
Goldreich & Scoville (1976) and of the OH maser emission from 
OH/IR stars by Elitzur, Goldreich, & Scoville (1976). The optically 
thick dusty shell surrounding the star was formed by a large mass 



loss rate of ~10~ 5 MQyr -1 . The shell expansion velocity of ~ 
20 km s _1 has been measured in OH maser and CO lines. Various 
methods and measurements suggest a distance to IRC+1001 1 in the 
range of 500 to 800 pc 

For an archetype of spherically symmetric AGB winds, the re- 
cent discovery by Hofmann et al. (2001; H01 hereafter) of distinct 
asymmetries in the IRC+10011 envelope came as a surprise. They 
obtained the first near infrared bispectrum speckle-interferometry 
observations of IRC+1001 1 in the J-, H- and K'-band with respec- 
tive resolutions of 48 mas, 56 mas and 73 mas. While the H- and 
K'-band images appear almost spherically symmetric, the J-band 
shows a clear asymmetry. Two structures can be identified: a com- 
pact elliptical core and a fainter fan-like structure. H01 also per- 
formed extensive one-dimensional radiative transfer modelling to 
explain the overall spectral energy distribution (SED) and angle- 
averaged visibility curves. Their model required a dust shell with 
optical depth r(0.55/i?n) = 30 around a 2250 K star, with a dust 
condensation temperature of 900 K. This one-dimensional model 
successfully captured the essence of the circumstellar dusty en- 
vironment of IRC+10011 but could not address the observed im- 
age asymmetry and its variation with wavelength. In addition, the 
model had difficulty explaining the far-IR flux, requiring an un- 
usual transition from a 1/r 2 density profile to the flatter 1/r for 
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r larger than 20.5 dust condensation radii. Finally, the model pro- 
duced scattered near-IR flux in excess of observations. 

We report here the results of 2D radiative transfer modelling 
of IRC+10011 that successfully explain the observed asymmetries. 
After analyzing in §2 general observational implications we de- 
scribe in §3 our model for a bipolar outflow in IRC+10011. In §4 
we present detailed comparison of the model results with the data 
and resolution of the problems encountered by the ID modelling. 
The discussion in §5 advances arguments for the role of bipolar 
jets in shaping the circumstellar envelope of IRC+1001 1 and other 
AGB stars. We conclude with a summary in §6. 



2 OBSERVATIONAL IMPLICATIONS 

The near-IR images, especially the J-band, place strong constraints 
on the dust density distribution in the inner regions. Emission at 
the shortest wavelengths comes from the hottest dust regions. For 
condensation temperature ~ 1,000 K the peak emission is at ~ 
4/xm, declining rapidly toward shorter wavelengths. At 1.24 fim, 
the J-band is dominated by dust scattering. It is easy to show that 
scattering by a l/r p dust density distribution produces a l/r p+1 
brightness profile (Vinkovic et al, 2003). The J-band image from 
H01 is elongated and axially symmetric. The decline of brightness 
from its central peak along this axis of symmetry is different in 
the opposite directions. In one direction it declines as 1/r 3 , corre- 
sponding to the 1/r 2 density profile typical of stellar winds. But 
in the other direction the brightness falls off only as 1/r 1 ' 5 , corre- 
sponding to the flat, unusual 1/r density profile. 

The large scale structure is not as well constrained by imaging. 
However, all observations are consistent with the following simple 
picture: An optically thick spherical wind has the standard 1/r 2 
density profile. Since the buildup of optical depth is concentrated in 
the innermost regions for this density law, the near-IR imaging pen- 
etrates close to the dust condensation region. The wind contains an 
imbedded bipolar structure of limited radial extent and density pro- 
file 1/r 0,5 . The system is observed at an inclination from the axis 
so that the wind obscures the receding part of the bipolar structure, 
creating the observed asymmetry of the scattering image, which 
traces directly the density distribution. The inclination angle must 
be < 45° since a larger value starts to expose the receding part. 
But the inclination cannot be too small because the approaching 
part would get in front of the wind hot dust, leading to a strong 
10 jj,m absorption feature, contrary to observations. Because of its 
shallow density profile, the column density of the bipolar struc- 
ture increases as r 5 away from the condensation cavity, and the 
size of J-band image corresponds to the distance where the scat- 
tering optical depth reaches unity. Regions further out do not show 
up because of self-absorption. Dust emission is affected also by the 
temperature distribution, and the central heating by the star tends to 
produce spherical isotherms. Images taken at longer wavelengths, 
such as the K-band, can thus appear more symmetric. 

Some qualitative estimates of the gas density follow imme- 
diately. The wind optical depth at the J-band must be > 1. This 
optical depth is accumulated close to the dust condensation radius, 
roughly 3 x 10 14 cm for a distance of 650 pc. Assuming a standard 
dust-to-gas mass ratio of 1 : 100, the gas density at the condensation 
radius is > 3xl0 7 cm -3 . For the bipolar structure, the J-band op- 
tical depth is ~ 1 across the size of the observed image, which is ~ 
2x 10 15 cm. This leads to a density estimate of ~ 7 x 10 6 cm~ 3 at 
the condensation radius within the bipolar structure. These rough 
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Figure 1. Sketch of the 2D model for the circumstellar dusty shell around 
IRC+10011. In a spherical wind with the standard 1/r 2 density profile are 
imbedded two polar cones with half-opening angle 6 C one and a 1/r - 5 
density profile. The system is viewed from angle i to the axis. 

estimates are within a factor 10 of the results of the detailed mod- 
elling described below. 

The density at the base of the outflow is about an order of 
magnitude lower in the bipolar structure than in the wind region. 
An outflow can bore its way through another denser one only if its 
velocity is higher so that it plows its way thanks to its ram pres- 
sure. The propagation of such high-velocity bipolar outflows has 
been studied extensively in many contexts, beginning with jets in 
extragalactic radio sources (Scheuer 1974). The jet terminates in a 
shock, resulting in an expanding, elongated cocoon similar to the 
observed bipolar structure. With a 1/r ' 5 density law, most of the 
bipolar structure mass is concentrated at its outer edge with the 
largest r, consistent with the structure of the expanding cocoon. 



3 2D MODELING OF IRC+10011 

Detailed 2D radiative transfer modeling was conducted with our 
newly developed code LELUYA (see Appendix). Starting with the 
minimal configuration that explains all available observations, we 
show on physical grounds that this minimal model must be modi- 
fied and propose a suitable modification. 

3.1 A Minimal Model 

For the first working model we adopt the geometry shown in figure 
Q which requires the minimal number of free parameters. Each po- 
lar cone is described by its half-opening angle # C onc and radial ex- 
tent -Rcono- Apart from discontinuities across the cone boundaries, 
the density depends only on r. It varies as 1/r ' 5 inside the cones 
and 1/r 2 outside, out to some final radius R ou t- To complete the 
description of the geometry we need to specify its inner boundary, 
and it is important to note that this cannot be done a-priori. Dust 
exists only where its temperature is below the condensation tem- 
perature T c . Following H01 we select T c = 900 K. The dust inner 
boundary, corresponding to the radial distance of dust condensa- 
tion, R c , is determined from 

T(R C (6))=T C (1) 

The equilibrium dust temperature, T, is set by balancing its emis- 
sion with the radiative heating. But the latter includes also the dif- 
fuse radiation, which is not known beforehand when the dust is 
optically thick; it can only be determined from the overall solu- 
tion. Furthermore, because the spherical symmetry is broken by the 
cones, the shape of the dust condensation surface can be expected 
to deviate from spherical and is not known a-priori. Therefore equa- 
tion Q completes the description of the geometry with an implicit 
definition of the inner boundary R c (8). 

The radiative transfer problem for radiatively heated dust pos- 
sesses general scaling properties (Ivezic & Elitzur 1997). As a re- 
sult, T c is the only dimensional quantity that need be specified. All 
other properties can be expressed in dimensionless terms. Luminos- 
ity is irrelevant, the only relevant property of the stellar radiation is 
its spectral shape, which we take as black-body at T* = 2,250 K. 
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Figure 2. SED observations and modeling. Data (see HOI) are indicated with various symbols and lines. Thick, smooth lines show the SED produced by two 
models: the dashed line corresponds to the minimal model ( i|3.11 . the solid line to the final, physical model ( i|3.21 . The inset shows an expanded view of the 
10pm region. All subsequent figures refer to the final physical model, whose properties are summarized in Table 1. 




Figure 3. Left: The model SED and its breakup to the stellar, dust scattering and emission components, as indicated. Right: Wavelength variation of the 
relative contribution of each component to total flux. Note the fast change from scattering to emission dominance around 2^tm. This transition is responsible 
for the observed wavelength variation of the image asymmetry in the near-IR. 



For individual dust grains, the only relevant properties are the spec- 
tral shapes of the absorption and scattering coefficients. For these 
we adopt spectral profiles corresponding to the silicate grains of 
Ossenkopf, Henning & Mathis (1992) with the standard size dis- 
tribution described by Mathis, Rumpl & Nordsieck (1977; MRN). 
Our calculations employ isotropic scattering. These properties are 
the same everywhere. 

Density and distance scales do not enter individually, only in- 
directly through overall optical depth. With two independent den- 
sity regions, the problem definition requires two independent op- 
tical depths. For this purpose we choose Tv and Ty, the overall 
optical depths at visual wavelengths along the axis and the equator, 



respectively. Spatial dimensions can be scaled with an arbitrary pre- 
defined distance, which we choose as the dust condensation radius 
in the equatorial plane, i? c (90°). Radial distance r is thus replaced 
with p = r/R c (90°) so that, e.g., p out = R oll t / Rc(90°). Equa- 
tionQbecomes an equation for the scaled boundary of the conden- 
sation cavity. The relation between angular displacement from the 
star $ and the distance p is 



2pZ 



(2) 



where is the stellar angular size and = i?*/i? c (90°) is the 
scaled stellar radius. Physical dimensions can be set if one specifies 



© 2004 RAS, MNRAS OOO.miTol 



4 Vinkovic et al. 



logpntenaity [Jy/orcsec 2 ]) log(lnten a ity [Jy/orcsec 2 ]) log(lntensity [Jy/orcsec 2 ]) 

8.5 9 9.5 9 9.5 10 9.5 10 10.5 



_ 




200 300 



rel. position (mas) rel. position (mas) rel. position (mas) 



Figure 5. Theoretical J-band (1.24/im), H-band (1.65/^m), and K'-band (2.12^tm) images of IRC+1001 1. Upper row: images for perfect resolution, without 
PSF convolution. The dot at the center of each image is the star. The nearby bright fan-shaped structure is scattered light escaping through the cone. Lower 
row: Images convolved with the instrumental PSF of H01 . Contours are plotted from 1 .5% to 29.5% of the peak brightness in steps of 1 %. The transition from 
scattered light dominance in the J-band to thermal dust emission in the K'-band creates a sudden disappearance of the image asymmetry. 



a stellar luminosity L, which determines the condensation radius 
i? c (90°). 

To summarize, in all of our model calculations the follow- 
ing quantities were held fixed: grain properties, T c = 900 K, T* 
= 2,250 K and the outer boundary p ou t = 1000. We varied Ty, 
i"V' 8 cone and pcone- Once a model is computed, comparison with 
observations introduces one more free parameter, the viewing an- 
gle i. A detailed discussion of the data is available in H01. Figure 
|2| shows the SED for the best fit model which has = 40, Tv = 
20, pcone = 700, and 8cone = 15°. The 10pm region is difficult 
to fit in full detail. Any further improvement would probably re- 
quire more complicated geometry and/or modified dust properties. 
The fit yields a bolometric flux of i-boi = 10 -9 W/m 2 , corre- 
sponding to = 10.82 mas for the stellar angular size, similar to 
the 10.9 mas derived in H01. The IR-fiux measurements determine 
the total amount of emitting dust. Assuming a standard ridCd/n = 



is 0.13 M 



0- 



3. 1. 1 Shortcomings of the Minimal Model 

The minimal-model success in fitting the SED implies that it con- 
tains the proper amount of dust. In addition to the SED, this model 



reproduces adequately all the imaging observations at near-IR. 
However, these observations probe only the innermost regions of 
the bipolar structure and do not provide adequate constraints on its 
full extent. 

We now show that the cones cannot extend all the way to p con e 
= 700, as required from the fit to the SED by the minimal model. If 
the cones were that large, the ratio of mass contained in them and 
in the wind region would be M cone /M w i n d = 1.7, that is, most of 
the circumstellar mass would be in the cones. But such large mass 
cannot be swept-up wind material because the fractional volume 
occupied by the cones is only 0.034. And building up this mass 
with enhanced outflow through the polar regions can also be ruled 
out, as follows: Mass conservation along stream lines yields v\t = 
R\ J Tjp 2 dp, where r)(p) = n(p)/ni is the dimensionless density 
profile, ui and n\ are the velocity and density at the streamline base 
Ri, and t is the duration of the outflow. Applying this relation to 



streamlines in the cone (rj = p 
yields 



-1/2 



) and wind (rj = p ) regions 



Vlt l _ 2 5/2 Vlt I 

D rPconc) r, 

t\\ I cone O Jtl I wind 



— Pout 



(3) 



With p C ono = 700, the product vit/Ri is 5.2xl0 6 in the cone re- 
gions while in the wind it is only 1000. Since the wind starts with 
a sonic velocity «i w ~ 1 km s - , the conical outflow would have 
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Figure 4. Visibility functions. Lines are model predictions, symbols are 
data points from HOI (near-IR) and Lipman et al. 2000 (11 pm). 



long as its temperature distribution corresponds to far-IR wave- 
lengths (equation A7 in Vinkovic et al. 2003). With the total far-IR 
flux as the only observational constrain, the only firm limit on this 
cold dust component is that it starts at p > 100 so that its tempera- 
ture is < 100 K; in all other respects, the geometry is arbitrary. 

To account for the far-IR excess, the H01 wind model, which 
did not incorporate the bipolar component, placed additional cold 
dust in a spherical component with p~ density profile at p > 20. 
Here we propose a simpler alternative: a detached spherical shell 
of increased dust density due to an earlier phase of higher mass 
loss rate. The radial wind dust density profile jumps by factor 3 at 
p = 100 so that 



?7wind 



1 for p < 100 

3 for 100 < p < 1000 



(4) 



With the cold dust displaced from the cones to the wind, the equa- 
torial optical depth increases from Ty = 40 in the minimal model 
to 40.72. Along the axis, r v drops from 20 to 5.1, of which 4.3 
comes from the cones. As is evident from figure |2| the SED pro- 
duced by this physical model is almost identical to that of the min- 
imal model. The contributions of different components to the total 
flux are shown in the left panel of figure|5| with the fractional con- 
tributions shown in the right panel. 

Table 1 summarizes the input parameters, and various proper- 
ties derived below, of our model. The proposal of a two-shell model 
is motivated mostly by physical plausibility since other than the 
SED, current observations do not place meaningful constraints on 
the cold dust configuration. We have verified that a disk geome- 
try for the cold dust would also successfully reproduce the SED 
by modeling with disk structures of p = 100 inner radius and ra- 
dial density profiles p -1 ^ 2 and p~ 2 . However, a disk geometry for 
the cold dust component suffers from the same shortcomings as the 
extended cones of the minimal model ( 83. Lit . 

The actual dust configuration is probably more complicated 
than our simple description. Imaging observations at 8.55 pm with 
spatial resolution of p ~ 50 by Marengo et al. (1999) suggest an 
extension along an axis almost perpendicular to the symmetry axis 
of the bipolar structure. This possible additional asymmetry is not 
accounted by our physical model. 



to start with velocity «i c ~ 5.2 x 10 3 ty,/t c km s _1 , where t„ and 
t c are the wind and cones lifetimes. This is impossible since the 
bipolar structure would extend much further than the wind even for 
t c = i w ; taking a physical t c <C t w only makes things worse. This 
argument can be easily extended to show that, irrespective of the 
magnitude of p con e, the mass in the cones could not be deposited 
purely by recent enhancement of polar mass loss rates. A substan- 
tial fraction, perhaps even all, of this mass must be swept-up wind 
material. 

3.2 A Physical Model 

While our minimal model provides successful fits for all observa- 
tions, pcono cannot be as large as this model requires. Unfortu- 
nately, observations do not yet meaningfully constrain p C one be- 
cause at p > 10 the near-IR brightness drops below current de- 
tection capabilities. A physical model for the origin of the bipolar 
structure ( 85 .21 suggests that the cones extend to p C onc = 47 (equa- 
tion[5}. Adopting this value, the minimal model must be modified 
to account for the far-IR flux produced by the large mass removed 
from the rest of the cones. This mass can be placed elsewhere as 



4 VISIBILITY FUNCTIONS AND IMAGES 

With the model parameters set from the SED, the surface bright- 
ness distribution is fully determined, and the visibility functions are 
calculated from the brightness. For comparison with observations, 
the visibility must be normalized with the flux collected within the 
field of view Gfv- If the image is divided into NxN pixels then 
the spatial frequency is qi = j/0fv, where i = 1...N. Adopting 
N — 300, sufficiently large to resolve the highest measured spa- 
tial frequencies, and Ofv from the instrumental system, the results 
shown in figure|4|contain no additional free parameters. In contrast 
with the SED, the visibility displays a strong sensitivity to the grain 
size. A change of only 0.05 pm in the maximum grain size a max 
has a significant effect on the visibility curves. Our physical model 
has a m ax = 0.20 pm, resulting in good fits for both the SED (figure 
|2j and the four different visibility curves (figure|4j. 

The J-band visibility is the most difficult to model because it is 
dominated by the scattered light and thus very sensitive to fine de- 
tails of the density distribution and grain size. Since the agreement 
between data and theory is better for small scales (higher spatial 
frequency), the quality of the fit to the J-band image can be ex- 
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stellar temperature 


2250 K 


stellar size 


10.8 mas (7 AUf) 


luminosity! 


1.3-10 4 L 


bolometric flux 


10 -9 Wm" 2 


dust condensation temperature 


900 K 


condensation radius i? c (90°) 


35.4 mas (23 AUf) 


viewing angle i 


25° 


Wind — Inner shell: 




density profile 


r -2 


inner boundary 


R c (6) (eq. 1) 


outer boundary 


100i? c (90°) 


radial Ty 


39.6 


mass* 


0.005 M Q 


ageft 


550 years 


mass loss rate* ft 


9-10~ 6 M Q yr- 1 


Wind — Outer shell: 




density profile 


r -2 


inner boundary 


100R C (90°) 


outer boundary 


1000R C (90°) 


radial ry 


1.1 


mass* 


0.13 Mq 


ageff 


5500 years 


mass loss rate* ft 


2-10~ 5 M Q yr" 1 


Cone Properties: 




opening angle 29 CO ne 


30° 


density profile 


r -0.5 


outer boundary 


47i? c (90°) 


axial tv 


4.3 


mass* 


10~ 4 Mq 


age*f 


< 200 years 



Table 1. Model parameters and derived properties. Quantities marked by f 
assume a distance of 650 pc, subject to an uncertainty of ±150 pc. Quanti- 
ties marked by tt additionally assume a wind velocity of 20 km s ~ 1 . Quan- 
tities marked by * assume n^a^/n = 10~ 21 cm 2 . Uncertainties and the 
acceptable range of the various properties are discussed in the text. 
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Figure 6. Temperature distribution around the condensation cavity. The 
contours start at 850 K and decrease at 50 K intervals. The dust condensa- 
tion temperature is 900 K. 

pected to deteriorate with distance from the star. The model does 
not explain the puzzling drop in the H-band visibility at q > 14 
cycles per arcsec, corresponding to structure smaller than the con- 
densation cavity. Since a similar drop is not present in the J-band, 
it must correspond to material that emits but does not scatter light 
significantly. Hot gas might be a possible explanation. 

Our model images and their convolution with the instrumen- 
tal PSF of H01 are shown in figure|5| The comparison between the 
model and observed images is satisfactory, indicating that the over- 
all geometry is properly captured by our simple model. The "halo" 
around the star in J-band model image is brighter than observed. 
Possible explanations are dust accumulation close to the equatorial 
region as well as asymmetric dust scattering. The overall image 
asymmetry is much more prominent in the J-band, where dust scat- 
tering dominates the radiative transfer (see figure As the wave- 
length shifts toward dominance of dust thermal emission, the image 
becomes more symmetric. The reason is that scattered light traces 



directly the density distribution while the dust emission is affected 
also by the temperature distribution. Figure [fj] shows the dust tem- 
perature distribution around the condensation cavity. Because of 
the central heating, the temperature decreases with radial distance 
and tends to create circularly symmetric isotherms. The asymmet- 
ric diffuse radiation distorts the isotherms, but the deviations from 
circularity are small, especially at the high dust temperatures traced 
by the K-band image. As a result, the image becomes more sym- 
metric, especially after convolution with the PSF as shown in the 
lower panel of figure|5| 

As evident from figure|5| the PSF convolution smears out the 
star and the nearby fan-shaped structure into one broad elongated 
peak whose center is shifted from the stellar position. This shift is 
more clearly noticeable in the brightness profiles, shown in figure 
Q The shift is 8.3 mas along the major axis in the J-band and 2.8 
mas for the H- and K'-bands. The images provide tight constraints 
on the inclination angle. Neither i = 20° nor i = 30° produce 
acceptable fits, so that i — 25° ± 3°. 



5 DISCUSSION 

Thanks to the scaling properties of dust radiative transfer, neither 
luminosity, distance or density absolute scales were specified. The 
distance to the source of 650±150 pc fixes those scales, so that 
the luminosity is 1.3 x 10 4 Lq and the dust condensation radius is 
R c (90°) = 23 ± 5 AU. The wind inner shell extends to 2,300 
AU, thus containing all the masers, including the OH 1612 MHz 
(Elitzur, Goldreich & Scoville 1976), and its mass is 5 x 10~ 3 M . 
With a wind velocity of 20 km s~ 1 , the duration of this phase is 550 
years and the corresponding mass loss rate is 9xl0~ 6 M yr -1 . 
The wind outer shell extends to 23,000 AU and its mass is 0.13 
Mq. Assuming the same wind velocity, the duration of this phase 
is 5,500 years and the corresponding mass loss rate is 2x 10~ 5 M 

yr- 1 . 

Modelling uncertainties allow a few times larger size of the in- 
ner shell, resulting in its larger mass and duration. The outer shell 
has much smaller overall uncertainties since its total mass is con- 
strained by the far-IR flux. We can still derive general conclusions 
that: 

(i) the duration of the inner shell is < 1,000 years, while the 
outer shell is 5 to 10 times older; 

(ii) the overall dust opacity comes mostly from the inner shell; 

(iii) the mass is contained almost exclusively in the outer shell, 
with about 10 to 30 times more mass than the inner shell; 

(iv) the mass loss rate is of the order of 10~ 5 Mq yr -1 , with a 
few times larger rate in the other shell; 

(v) the overall circumstellar mass of 0.13 Mq indicates that 
IRC+1001 1 is close to the end of its AGB evolution. 



5.1 Dust Properties 

Our models employ silicate grains from Ossenkopf et al. (1992) 
with the standard MRN size distribution. We found that the up- 
per limit on the grain sizes had to be reduced to a max = 0.20 /im 
from the standard 0.25 ^m. While this change made little differ- 
ence in the SED analysis, it was necessary for proper fits of the vis- 
ibility curves. The most important effect of a max is control of the 
crossover from scattering to emission dominance, crucial for expla- 
nation of the observed change from elongated to circular images 
between the J- and K-bands (see f0J. Although we cannot claim to 
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Figure 7. Brightness profiles along the major and minor axes in J (top 
two panels), H (middle two panels) and K (bottom two panels). Thick lines 
show the model predictions with and without PSF convolution. The thin 
lines show the profiles from the H01 data above the noise level (within 1.5% 
of the peak brightness). The strong central peak in the theoretical profile is 
the star, while the secondary peak visible on the major axis is light scattered 
from the polar cone. 



have determined the precise magnitude of a max , the fact that it is 
smaller than the standard seems certain. 

The dust properties in our model were the same everywhere 
to minimize the number of free parameters. In a detailed study of 
the proto-planetary nebula IRAS 16342-3814, Dijkstra et al. (2003) 
find that the maximum grain size varies from ~ 1.3 in a torus 
around the star to ~ 0.09 ^tm in the bipolar lobes. If such a variation 
in dust properties can occur already on the AGB, the a max we find 
would represent an average over the cones and wind regions. 

5.2 Jet Model for the Bipolar Structure 

The near-IR brightness observations map a region of the cones that 
extends to p ~ 8 and has optical depth -ry ~ 1.4, corresponding 
to a gas density at the base of each cone of m c = 1.3 x 10 6 cm -3 . 
In contrast, the gas density at the base of the wind region (obtained 
from Ty = 41) is ni w = 1.7 x 10 s cm -3 . The large density disparity 
amplifies our earlier conclusion that the bipolar cones are sustained 
by high-velocity ram pressure. 

The small density at the base of the cones shows that their ma- 
terial has been evacuated and deposited at larger distances by a re- 
cent event. We propose the following simple scenario for the bipo- 
lar structure: High-velocity low-density jets were recently turned 
on at the polar regions. The jets cleared out polar cavities but are 
trapped by the material pushed ahead by their ram pressure, result- 
ing in an expanding cocoon as described first by Scheuer (1974). 
Our model cones are a description of the current density distribu- 
tion of the cocoon, a snapshot of an inherently dynamic structure. 
In this picture, the mass in the cones is swept-up ambient wind ma- 
terial and the cone boundary is then 

/ E \ 2 / 3 

P cone=(^) =47. (5) 

\2 nic / 

The swept-up mass is only ~ 10~ 4 M@. The leading edge of the 
cocoon moves at velocity v c = /3v„, where v w is the local wind 
velocity and /3 > 1. From pressure balance during jet confinement, 

n c (v c - v„) 2 = n w Vt (6) 

where n c and n w are densities across the cocoon leading edge and 
v t is the local speed of sound in the wind. This condition requires 
that the density of the cones be smaller than the ambient density 
into which they are expanding, i.e., n c < n w , restricting the cone 
radial extension to p 26 which is slightly smaller than the de- 
rived pcone- We attribute this discrepancy to the approximate nature 
of our model, in which the complex structure of the cocoon-wind 
boundary is replaced with the sharp-cutoff of the simple power- 
law density distribution of the cones. Taking n c ~ n w at the cone 
boundary, pressure balance implies v c ~ w w , consistent with a re- 
cent start of the jet confinement. Assuming that the cocoon radial 
boundary moves according to p CO nc cc t a , with a > 1 to ensure ac- 
celeration, its velocity is v c — ap conc Ri jt. This yields an estimate 
for the jet lifetime 

*jet = % — Pconc < 200 years (7) 

p «w 

fora/P < 1. 

Because of the steep decline of the wind density, the expan- 
sion accelerates rapidly as the cocoon boundary reaches lower den- 
sity regions. Eventually it will break out of the wind, exposing the 
underlying jets. Indeed, a striking example of such a configuration 
comes from the recent observations, including proper motion mea- 
surements, of water masers in W43A by Imai et al. (2002). The 
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Figure 8. Schematic representation of our proposed bipolar jet evolution 
during the late stages of AGB winds. The top panel shows an early phase 
when the jets are still confined by the spherical wind. IRC+1001 1 is the ear- 
liest manifestation of this structure. The bottom panel illustrates evolved jets 
that have broken out of the spherical confinement. CIT6 and W43A show 
clear evidence of such structure. Several other sources give indications of 
intermediate stages and are listed according to approximate ordering sug- 
gested by current observations. The source RV Boo (Bergman et al. 2000, 
Biller et al. 2003) gives an indication of bipolarity, which is less certain than 
in the other objects. 

observations reveal tightly collimated velocities of ~ 150 km s 
at distances up to ~ 0.3 pc at the two ends of an axis through 
the star. These masers are created by the impact of the jets on 
clumps in the surrounding medium. In addition to these far-away 
high-velocity masers, the source displays the usual configuration 
typical of OH/IR stars - OH and water masers in shells expand- 
ing with velocities ~ 9 km s _1 with radii of ~ 500 AU. Therefore 
this source displays both the spherical AGB wind and the jets that 
broke through it. A similar, and probably more evolved, example is 
IRAS16342-3814. Similar to W43A, this is a "water-fountain" jet 
PPN source but its bipolar cavities are also seen with HST (Sahai 
et al. 1999). 

5.3 Asymmetry Evolution in AGB Stars 

IRC+1001 1 and W43A can be considered, respectively, the 
youngest and most evolved examples of sources displaying the evo- 
lution of bipolar jets working their way through AGB winds. The 
proposal that jets, operating at the late AGB or early post- AGB 
phase, are the primary mechanisms for shaping PNe has been made 
already by Sahai & Trauger (1998) and since supported by numer- 
ous observations. The prototype C-rich star IRC+10216 shows cir- 
cular shape on the 20 arcsec scale both in V-band (de Laverney 
2003) and molecular line images (e.g., Dayal & Bieging 1995). 
But high-resolution IR imaging at the 0. 1 arcsec scale reveal elon- 
gated structure similar to that in IRC+1001 1 (Osterbart et al. 2000, 



Weigelt et al. 2002). Unlike IRC+1001 1, though, where only the J- 
band image gives clear indication of asymmetry, in IRC+10216 it is 
evident even in the K-band. This strongly suggests that IRC+10216 
represents a more advanced stage than IRC+1001 1 of the evolution 
of a jet-driven cocoon confined by the ambient spherical wind. 

The C-rich star V Hya provides an example that is further 
along in evolution. Recent CO observations by Sahai et al. (2003a) 
show that the bulk of the emission comes from an elongated struc- 
ture centered on the star. In addition, an emission blob is approach- 
ing at a projected line-of-sight velocity of 250 km s~ 1 along an axis 
perpendicular to this elongation. This is the expected morphology 
of a bipolar outflow breaking from the confinement of the high- 
density region of the AGB wind if the receding blob is obscured 
by the central torus. A similar structure has been found in the O- 
rich star X Her. Partially resolved CO observations by Kahane & 
Jura 1996 reveal a spherical component expanding with only 2.5 
km s _1 and two symmetrically displaced 10 km s _1 components, 
likely to be the red and blue shifted cones of a weakly collimated 
bipolar flow. The bipolar lobes are ~ 1.5 times bigger than the 
spherical component. Finally, the C-rich star CIT6 presents an even 
more evolved system. A bipolar asymmetry dominates the image in 
molecular line mapping by Lindqvist et al. (2000), Keck imaging 
by Monnier et al. (2000) and HST-NICMOS imaging by Schmidt 
et al. (2002). 

Figure|8|shows a schematic sketch of our proposed evolution- 
ary sequence. The evolutionary stage of each indicated object is 
our rough estimate based on current observations. This figure is 
only meant as an illustration of the time-line suggested by our pro- 
posed scenario. The placing of different objects is based on differ- 
ent kinds of data. For example, for X Her we only have single-dish 
mm-wave observations with their attendant large beams, whereas 
IRC+10216 and V Hya were observed with the much higher reso- 
lution of HST. Also, whereas there is direct kinematic evidence for 
the high-velocity jets in V Hya and W43A, the same is lacking for 
the other objects. We can also expect that in objects with different 
jet properties (e.g., mass flux, speed, opening angle) and different 
AGB mass-loss rates in the inner region, the jets will follow differ- 
ent time histories of when they "break out" and the opening angle 
of the bipolar cone which they dig in the AGB envelope will be 
different. Also, if jets are episodic than they can change their direc- 
tion (for which there is observational evidence). The actual picture 
can be expected to be quite more complex than the simple sketch in 
Fig [HI Nevertheless, we expect the displayed sequence to provide 
useful guidance for future studies. 



6 SUMMARY AND CONCLUSIONS 

We find that the circumstellar shell of IRC+1001 1 contains ~ 0.13 
M , extends to a radial distance of ~ 23,000 AU (~ 35") from 
the star and is ~ 5,500 years old. Most of the mass (~ 96%) is 
contained in the outer shell from ~ 2,300 AU (~ 3.5"), corre- 
sponding to an earlier phase when the mass-loss rate was about 
factor 3 higher than now. The near-IR image asymmetries discov- 
ered within the central ~ 0. l"of this system originate from ~ 10~ 4 
Mq of swept-up wind material in a cocoon elongated along the 
axis, extending to a radial distance of ~ 1,100 AU. The cocoon 
expansion is driven by bipolar jets that it confines, and that were 
switched on < 200 years ago. The axial symmetry of the J-band 
image eliminates the possibility of a companion star, unless closer 
than ~ 5 stellar radii. Higher sensitivity and/or better angular reso- 
lution would uncover image asymmetry in the K-band too. 
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Jet-driven cocoon expansion at various stages of development 
has now been observed in a number of AGB stars, culminating in 
breakout from the confining spherical wind ( £I5.3> . The immedi- 
ate post-AGB stage is believed to be the proto-planetary-nebula 
(PPN) phase. A morphological study of a large sample of PPNs 
suggested the presence of jets that broke through the massive AGB 
wind, but it was not known if a similar morphology extends back to 
the AGB phase (Meixner et al. 1999; Ueta, Meixner & Bobrowsky 
2000). Indeed, jets are found to be quite common in PPN as shown 
by the recent observations of IRAS 16342-3814 (Sahai et al. 1999) 
K3-35 (Miranda et al. 2001), Hen 3-1475 (Riera et al. 2003) and 
IRAS22036+5306 (Sahai et al 2003b), for example. The case of 
K3-35 is particularly striking because of its great similarity to the 
AGB star W43A: water masers at the tips of bipolar jets at a large 
distance from the systemic center, which is surrounded by masers 
in the standard spherical shell configuration. This strongly suggests 
that W43A provides a glimpse of the immediate precursor of K3- 
35. 

These new developments enable us to identify the first in- 
stance of symmetry breaking in the evolution from AGB to plan- 
etary nebula. Bipolar asymmetry appears during the final stages of 
AGB mass outflow. Mounting evidence suggests that this asymme- 
try is driven by collimated outflow in the polar regions. More com- 
plex geometries emerge in the post-AGB phase from a mixture of 
various processes that could involve multiple jets, fast winds, etc. 
These processes operate in the environment shaped by the AGB 
phase, leading to the myriad of complex structures found in PPN 
sources (e.g. Su, Hrivnak & Kwok 2001). 
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APPENDIX A: DETAILS OF RADIATIVE TRANSFER 
MODELING 

LELUYA (www.leluya.org) is our newly developed 2D radiative 
transfer code that works with axially symmetric dust configura- 
tions. It solves the integral equation of the formal solution of ra- 
diative transfer including dust scattering, absorption and thermal 
emission. The solution is based on a long-characteristics approach 
to the direct method of solving the matrix version of the integral 
equation (Kurucz 1969). The equations are solved on a highly un- 
structured triangular self-adaptive grid that enables LELUYA to re- 
solve simultaneously many orders of magnitude in both spatial and 
optical depth space. It also enables automatic reshaping of the dust- 
free cavity around the central source according to asymmetries in 
the diffuse radiation (eq. 1). All grid points are coupled with each 
other through a correlation matrix based on the dust scattering. A 
simple matrix inversion determines the solution of radiative trans- 
fer for a given dust temperature distribution without any iterations. 
The temperature is then updated and the procedure repeated. Lumi- 
nosity conservation within 5% is achieved in only three steps. 

Figure \Al\ shows LELUYA's computational grid for the best 
fit minimal model described in A3, II Three grids of different res- 
olutions were created for three sets of wavelengths, based on the 
density and optical depth variation. The first grid has 2982 points 
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Al Luminosity Conservation 

Luminosity conservation is the test determining convergence to the 
correct physical solution. A decrease in computed luminosity indi- 
cates energy sink due to insufficient spatial grid resolution, while an 
increase reflects energy excess due to a coarse angular grid. It is im- 
portant to note that because of the lack of spherical symmetry, the 
bolometric flux does vary over spherical surfaces. The conserved 
quantity is luminosity, the energy transmitted per unit time across 
any surface enclosing the star. For a sphere of radius p, the lumi- 
nosity is computed from the radial component of the bolometric 
flux vector Fboi r via 



L(p) = Anp 2 [ F hol>I (p,9)d cos 9 
Jo 

and the luminosity conservation relation is 

Hp) _ , 



(Al) 



(A2) 



at every p, where L is the stellar luminosity. In spherical sym- 
metry Fboi, r is ^-independent and 47rp 2 _Fboi,r/L = 1. When 
the spherical symmetry is broken Fb a \ t i becomes 9-dependent 
and 4Tvp 2 F hol , T (9)/h can exceed unity in certain directions, cor- 
responding to locally enhanced energy outflow. 

Our model calculations conserve luminosity within 5% at all 
radii. Figure lA2l shows the angular variation of Fboi,r{9) and its 
following five contributions: stellar, inward and outward emission, 
and inward and outward scattered flux. These angular variations are 
shown at p = 1.1, 500 and 1000. The small spikes in .Fboi.r close 
to Scone are real, reflecting the irregular shape of the dust conden- 
sation surface. Even though these irregularities are spatially small, 
their effect on optical depth variations magnifies their importance. 
At small radii, energy outflow through the cones is enhanced in 
comparison with the wind and is the main reason for their higher 
temperature. This region is dominated by stellar contribution. At 
large radii these roles are reversed, the diffuse radiation (mostly 
dust emission) takes over and the temperatures inside and outside 
the cones become equal. 



Figure A2. Angular dependence of the radial bolometric flux over spheres 
of radius p = 1.1, 500 and 1000 for the best fit minimal model described in 
CHI The numerical precision of luminosity conservation (eq, IA2l is indi- 
cated from the listed L(p) /L in each panel. 



and starts with r c = 120 at 0.2pm, the shortest wavelength con- 
sidered; this is the grid shown in the figure. The second grid has 
2836 points and starts at wavelengths with r c = 1.2. The third has 
2177 grid points for wavelengths with r c ^ 0.1. Angular integra- 
tion around a grid point is performed over a highly non-uniform 
self-adaptive angular grid (with about 550 rays on average). 
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